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Abstract 

In this paper an algorithm is given to determine all possible structurally different 
linearly conjugate realizations of a given kinetic polynomial system. The solution is 
based on the iterative search for constrained dense realizations using linear program¬ 
ming. Since there might exist exponentially many different reaction graph structures, 
we cannot expect to have a polynomial-time algorithm, but we can organize the com¬ 
putation in such a way that polynomial time is elapsed between displaying any two 
consecutive realizations. The correctness of the algorithm is proved, and possibilities 
of a parallel implementation are discussed. The operation of the method is shown 
on two illustrative examples. 

Keywords; reaction networks, reaction graphs, linear conjugacy, linear programming 


1 Introduction 

Chemical reaction networks (CRNs) obeying the mass action law can be originated from 
the dynamical modelling of chemical and biochemical processes. However, these models 
are also capable of describing all important phenomena of nonlinear dynamical behaviour 
and have several applications in different fields of science and engineering OEl. Kinetic 
systems, which describe the dynamical behaviour of the CRNs have a simple algebraic 
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characterization, which makes it possible to develop effective computational methods for 
dynamical analysis, and even control mmu- 

The graph representation of reaction networks, especially the exploration of the rela¬ 
tion between the reaction graph structure and the dynamics of the network (preferably 
without the precise knowledge of the reaction rate coefficients) has become an impor¬ 
tant research area in chemical reaction network theory since the 1970’s. One of the first 
comprehensive overviews of mass-action-type CRNs can be found in clearly indicating 
that the class of such kinetic systems is far more general than the description of “chemi¬ 
cally reacting mixtures in closed vessels”. In the same article, the Kirchhoff matrix based 
description of reaction networks is also introduced. 

Complex balance is a fundamentally important property of CRNs (see, e.g miin]). 
Roughly speaking, complex balance means that the signed sum of the incoming and 
outgoing fluxes corresponding to each vertex (complex) at equilibrium is zero. It is an 
interesting and important property that the set of equilibrium points of any complex bal¬ 
anced system forms a toric variety in the state space nni. One of the most significant 
achievements in chemical reaction network theory is the Deficiency Zero Theorem saying 
that weakly reversible reaction networks of deficiency zero are complex balanced, inde¬ 
pendently of the values of the reaction rate coefficients m- An important recent result 
related to complex balance is the possible general proof of the Global Attractor Conjec¬ 
ture [12]. According to this conjecture, complex balanced CRNs with mass action kinetics 
are globally stable within the positive orthant with a known logarithmic Lyapunov func¬ 
tion that is independent of the reaction rate coefficients. The experimentally observed 
robust dynamical behavior of certain reaction mechanisms can also be understood from 
the special properties of the network structure |13j . 

It is known from the “fundamental dogma of chemical kinetics” that the reaction 
graph corresponding to a given kinetic dynamics is generally non-unique |m [T5] . This 
phenomenon is also called dynamical equivalence or macro-equivalence |6],[l6], and 
it has been studied in the case of Michaelis-Menten kinetics, too m- Clearly, this kind 
of structural non-uniqueness may seriously hamper the parameter estimation (inference) 
of reaction networks from measurement data, that is generally a challenging task [181119] . 
To improve the solvability of such inference problems, it is worth involving as much prior 
information as possible into the computation problem in the form of additional con¬ 
straints [201 EH [22]. Necessary and sufficient conditions for a general polynomial vector 
field to be kinetic were first given in [I6]. In the constructive proof, a procedure was 
given to construct one possible dynamically equivalent reaction network structure (called 
the canonical structure) realizing a kinetic dynamics. The first theoretical results about 
dynamical equivalence of CRNs were published in [23] pointing out the convex cone struc¬ 
ture of the parameter space. Motivated by these results, the numerical computation of 
dynamically equivalent realizations was first put into an optimization framework in m 
by proposing a method for determining the graph structures containing the minimum or 
maximum number of reactions. It is also important to know that several fundamental 
properties of CRNs such as (weak) reversibility, complex and detailed balance or defi¬ 
ciency are not ‘encoded’ into the kinetic differential equations but may vary with the 
dynamically equivalent structures corresponding to the same dynamics. This fact is an 
important motivating factor to develop methods for the computation of reaction graphs 
of kinetic systems. 

There is a significant extension of dynamical equivalence called linear conjugacy, where 
the kinetic model is subject to a positive definite linear diagonal state transformation |25) . 
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Obviously, linear conjugacy preserves the main qualitative dynamical properties of CRNs 
like stability and multiplicities or the boundedness of solutions. However, due to the larger 
degree of freedom introduced by the transformation parameters, it allows a wider variety 
of possible reaction graph structures than dynamical equivalence does. Therefore, several 
optimization based computational methods have been suggested to hnd linearly conjugate, 
or as a special case, dynamically equivalent realizations of kinetic systems having preferred 
properties such as density/sparsity, complex or detailed balance, or minimum dehciency 
|26) . An algorithm for hnding all possible sparse dynamically equivalent reaction network 
structures was proposed and applied for a Lorenz system transformed into kinetic form 
in |27], using the mixed integer linear programming (MILP) framework. 

After treating the above mentioned important special cases, a question arises naturally: 
Is it possible to give a computationally efficient algorithm for determining all possible re¬ 
action graph structures corresponding to linearly conjugate CRN realizations of a given 
kinetic system? The aim of this paper is to give a solution to this problem. Of course, we 
cannot expect to hnd a polynomial-time algorithm for the overall problem, since exponen¬ 
tially many different realizing graph structures may exist for a kinetic model. However, 
as it will be shown, it is possible to achieve that each computation step (i.e. hnding and 
displaying the next realization after the previous one) is performed in polynomial time. 

2 Basic notions 

In this section, we summarize the basic notions and results for both algebraic and graph- 
based representations of CRNs. We will use the following notations: 

M the set of real numbers 
M+ the set of nonnegative real numbers 
N the set of natural numbers, including 0 

j^nxm gg|- matrices having entries from a set H in n rows and m columns 
[M]ij the entry in row i and column j of matrix M 

2.1 Algebraic characterization 

It is important to remark here that similarly to im. n or j2], we consider reaction net¬ 
works as a general system class representing nonlinear dynamical systems with nonnega¬ 
tive states, therefore we do not require that they fulhl actual physico- chemical constraints 
like mass conservation. In other words, there is no constraint on how diherent complexes 
can transform into each other in the network. 

Definition 2.1. Chemical reaction networks can he determined by the following three sets 
(see e.g. [28l E]). 

• A set of species: S = {Xi | z G {1,...,«.}} 

• A set of complexes: C = {Cj | j G {1,..., m}}, where 

n 

= J'e {!,...,m} 

i=l 

aji G N j G {1,..., m}, z G {1,..., n} 

The complexes are formal linear combinations of the species with coefficients, called 

the stoichiometric coefficients. 
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• A set of reactions: IZ C {{Ci, Cj) \ Ci, Cj G C} 

The ordered pair {Ci,Cj) corresponds to the reaction Ci —)■ Cj. 

To each ordered pair {Ci, Cj) where i,j G {1,... m} and i ^ j there belongs a nonnegative 
real number kij called the reaction rate coefficient. The reaction Ci —)■ Cj takes place 
if and only if the corresponding coefficient is positive. 


The properties of the reaction network are encoded by special matrices. 

Definition 2.2. Y G is the complex composition matrix of the CRN if its 

entries are the stoichiometric coefficients. 

[Y]ij = Oji i e {l,...,n}, j e {l,...,m} (1) 

Definition 2.3. G is the Kirchhoff matrix of the CRN if its entries are 

determined by the reaction rate coefficients as follows: 




kji ifi^j 

m 

- kii if i= j 

1=1,i^i 


( 2 ) 


Since the sum of entries in each column is zero, Ak is also called a column conservation 
matrix. 


Let X : M —)■ Wf be a function, which defines the concentrations of the species 
depending on time. Assuming mass-action kinetics, the dynamics of the function x can 
be described by dynamical equations of the following form; 


x = Y ■ Ak- f){x) (3) 

where : Wf —)■ is a monomial-type vector-mapping, 

n 

= n j = (4) 

i=l 


According to Definitions 2.2, 2.3 and Equations (§. 0 the dynamics of a reaction 
network can be given by a set of polynomial ODEs, but not every polynomial system 
describes a CRN. 


Definition 2.4. Let x : M —)■ M" be a function, M G a matrix and ip : —)■ a 

monomial function. The polynomial system 


X = M ■ ip{x) (5) 

is called kinetic if there exist a matrix Y G and a Kirchhoff matrix Ak G 

so that 

M ■ ip{x) = Y ■ Ak - 'f{x) (6) 

where : M” —)■ is a monomial function determined by the entries of matrix Y, 

= nr=i j e {1,..., m}. 
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The matrices Y and Ak characterize the dynamics of the kinetic system, as well as the 
CRN. However, the polynomial kinetic system does not nniquely determine the matrices. 
Reaction networks with different sets of complexes and reactions can be governed by the 
same dynamics, see e.g. |6l |23l [23]. If the matrices Y and of a reaction network fnlfil 
Equation ([^, then the CRN is called a dynamically equivalent realization of the 
kinetic system ([^, and it is denoted by the matrix pair (R, Ak). 

The notion of dynamical equivalence can be extended to the case when the polynomial 
system is subject to a positive linear diagonal state transformation. It is known from j^H] 
that such a transformation preserves the kinetic property of the system. 

Let T G be a positive definite diagonal matrix. The state transformation is 

performed as follows: 

X = T ■ X, X = T~^ ■ X (7) 

Applying it to the polynomial system ([^ we get 

^ = T-^ ■ X = T-A M ■ (f(x) =T-^ -M ■ ip{T • x) = ■ M ■ • ip{x) (8) 

where G is a positive definite diagonal matrix so that = 7>i{T ■ 1 ) for 

i G {1,..., n}, is the ith coordinate function of 99 , and 1 G M” is a column vector with 
all coordinates equal to 1. Now we can give the extended definition. 

Definition 2.5. A reaction network (F, Aj,) is a linearly conjugate realization of the 
kinetic system ([^ if there exists a positive definite diagonal matrix T G so that 

Y -A^- ij{x) =T-^ -M -^T- ^{x) (9) 

where Y G ; MJf —)■ with ipjfx) = nr=i for j G {!,..., m}, and 

A'f, G is a Kirchhoff matrix. 


It can be seen that dynamical equivalence is a special case of linear conjugacy, when 
the matrix T, and therefore the matrices T~^ and 'h-r as well are identity matrices. 

Since the monomial functions cp and xj in Equations ([^ and (|^ might be different, in 
both cases the set of complexes is not fixed. By applying the method described in |16j a 
suitable set of complexes can be determined, but there are other possible sets as well. It 
is clear that the complexes determined by the monomials of function tp must be in the 
set C, but arbitrary further complexes might be involved as well, which appear in the 
original kinetic equations with zero coefficients. These additional complexes change the 
dimensions of the matrices Y and Aj., therefore we have to modify the matrices M and 
as well, in order to get the following equation; 


F ■ A'fc • x/j{x) = T-^ ■ M' ■ . ^jJ{x) 


( 10 ) 


where the matrices M' G and G 


have the same columns and diagonal 


entries as M and belonging to the complexes determined by (^, and zero columns and 
1 diagonal entries belonging to all additional complexes, respectively. 

Since there is the same monomial-type vector-mapping xj on both sides of Equation 
(10), the equation can be fulfilled if and only if the coefficients belonging the same mono¬ 
mials are pairwise identical. This means that by using the notation Ak = A'^ ■ we 

can rewrite Equation (10) as 
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Y ■Ak = T-^ ■ M' (11) 

where Ak is a Kirchhoff matrix, too, obtained by scaling the columns of A'j^ by positive 
constants. It is easy to see that this operation preserves the set of reactions, but changes 
the values of the non-zero entries. The actual reaction rate coefficients of the linearly 
conjugate network are contained in the matrix 

A't = At ■ 4'j. (12) 

From now on we will consider only linearly conjugate realizations on a hxed set of 
complexes. A reaction network which is a linearly conjugate realization of a kinetic system 
can be identihed by its matrices T, Y, and A'/^. However, since matrix Y is hxed, and the 
matrix Ak is returned by the computation, we will simply denote this realization by the 
matrix pair (T, Ak). 

2.2 Graph representation 

A reaction network can be described by a weighted directed graph. 

Definition 2.6. The graph G{V, E) representing the CRN is called Feinberg-Horn- 
Jackson graph, or reaction graph for short. If the weights are given by the function 
w : E{G) —)■ M+, the reaction graph is defined as follows: 

• the vertices correspond to the complexes, V(G) = C, 

• the directed edges describe the reactions, E{G) = IZ, 

there is a directed edge from vertex Gi to vertex Gj if and only if the reaction Gi —)■ Gj 
takes place (i.e. kij > 0), 

• the weights of the edges are the reaction rate coefficients, 
w{{Gi,Gj)) = kij where {Gi,Gj) G IZ. 

Loops and multiple edges are not allowed in a reaction graph. 

The applicability of the algorithm presented in this paper depends on a recently proved 
important property of the so-called dense realizations. Therefore, we formally define the 
notion of dense realizations, and then recall the related result from |3n) . 

Definition 2.7. A realization of a CRN is a dense realization if the maximum number 
of reactions take place. 

It is easy to see that the dense property of a CRN realization is equivalent to the feature 
that its Kirchhoff matrix Ak has the maximum number of positive off-diagonal entries. 
In a set of weighted directed graphs we call a graph super-structure if it contains each 
element as a subgraph not considering edge weights, and it is minimal under inclusion. It 
is clear that the structure of the super-structure graph is unique, because there cannot be 
two different graphs that are subgraphs of each other, but have different structures. The 
following proposition published in [20] says that dense linearly conjugate realizations have 
a super-structure property even if an arbitrary additional finite set of linear constraints 
on the rate coefficients and the transformation parameters has to be fulfilled. 
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Proposition 2.8. (SO] Among all the realizations linearly conjugate to a given kinetic 
system and fulfilling a finite set of additional linear constraints there is a realization 
determining a super-structure. 


An important example of the linear constraints mentioned in Proposition 2.8| is the 
condition for mass conservation. According to this, the total mass in the CRN is preserved 
(i.e. the system is called kinetically mass conserving [31]) if and only if there exists a 
strictly positive vector k G MA such that the following holds; 


k' - Y -A. = 0 


T 


(13) 


where 0 G M"’ is the null vector. The chemical meaning of the vector k is that its entries are 


the molecular (or atomic) weights of the species in the network. By using Equation (12) 
it is easy to see that for a given vector fc, Equation (13) holds if and only if {k~^ ■ Y ■ A'jf) 


is the null vector, since is an invertible diagonal matrix. Therefore, the computation 
of all reaction graph structures that describe realizations obeying mass conservation with 
a given vector k is also possible by using the method proposed in Section and adding 
Equation (13) to the computation constraints. 


3 Computation model for computing linearly conjugate 
realizations 

Linearly conjugate realizations can be computed by using a linear optimization model 

Equation 0 must be fulfilled. We assume that 


As it was described in Section 2.1 


the set of complexes is given, so the coefficient matrix does not need any modification, 
therefore M = M' holds. Consequently the equation characterizing linearly conjugate 
realizations can be written as follows: 


T-Y M -Y ■ A. = 0 


(14) 


where 0 G denotes the zero matrix. The matrices Y and M are fixed, and the 

variables are contained in the matrices T~^ and A^. Specifically, the variables are the 
off-diagonal entries of the matrix Ak and the diagonal entries of the matrix T~^. We do 
not consider the diagonal entries of matrix A^ as variables, because these are uniquely 
determined by the off-diagonal entries of matrix 




(15) 


i=i 

j¥=i 


Equation (14) guarantees linear conjugacy, and Equations (15), (16), and 0 ensure that 
the matrices T~^ and Ak meet their definitions. 


[Ak]ij>0 i,j e iy^ j (16) 

[T~%>0 iG {!,...,n} (17) 

In our algorithm it will be necessary to exclude some set "H C 77 of reactions from the 
computed realizations. This can be written in the form of a linear constraint as follows: 

[Akfi. = 0 iQ,Cfien (18) 
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The feasibility of the linear constraints (14)-(18) can be verihed, and valid solntions (if 
exist) can be determined in the framework of linear programming. 

It is important to remark that in the above compntation model, the bonndedness prop¬ 
erty of all variables can be ensnred so that the set of possible reaction graph strnctnres 
remains the same as it was proved in 


Proposition 3.1. [SOj For any linearly conjugate realization {T,Ak) of a kinetic system 
there is another linearly conjugate realization {T*,Al) with all variables smaller than 
the given upper bound(s) so that the reaction graphs describing the two realizations are 
structurally identical. 


In onr algorithm it is necessary to compnte constrained dense linearly conjngate re¬ 
alizations, which rises two technical problems. The hrst is that we have to maximize the 
nnmber of positive off-diagonal entries of matrix Ak- The second one is that there are 
strict ineqnalities (the diagonal entries of matrix T~^ mnst be strictly positive). There 
are several possibilities for handling these issnes. In onr implementation, we will apply 
the method presented in [30], which can determine constrained dense linearly conjngate 
realizations in polynomial time. The reasons for this choice are briefly the following. 
Several solntions nse mixed integer linear programming (MILP) for determining dense 
realizations (see e.g. [26] and the references therein). However, MILP problems are known 
to be NP-hard and their application wonld not allow the time complexity of the compn¬ 
tation between displaying any two consecntive realizations to be polynomial. Moreover, 
the treatment of strict ineqnalities in optimization is not trivial and needs special atten¬ 
tion [32]. Therefore, ineqnalities like x > c are often transformed in practice to the form 
X > c + e, where £ > 0 is a small nnmber, bnt this may make the compntations less 
accnrate. The method in (30) avoids both issnes mentioned above by ntilizing the spe¬ 
cial properties of linearly conjngate realizations, and therefore, this is the most reliable 
polynomial-time method that we cnrrently know. 


4 Algorithm for computing all realizations 


Linearly conjngate realizations are parametrically not nniqne, the matrices A^ and T~^ in 
Eqnation (14) can be scaled by the same arbitrary positive scalar [23113U] . Therefore, onr 
aim is to determine all the possible reaction graph strnctnres describing linearly conjngate 
realizations of a kinetic system. In this section we present and analyze an algorithm for 
compnting all snch reaction graph strnctnres. This is the main resnlt of the paper. 

From now on the reaction graphs will be assnmed to be nnweighted directed graphs, 
i.e. by reaction graph we will mean only its strnctnre. According to Proposition |2.8 


if Gd is the reaction graph of the dense realization and graph Gr describes another 
realization, then E{Gr) C E{Gd) holds. As it was proven in [ 30 ], the dense realization 
can be determined by a polynomial algorithm, which is the hrst step of onr method. 

There might be reactions that take place in each realization. These are called core 
reactions, and the edges representing them are the core edges. The set of core edges 
is denoted by Ec, which can also be determined by a polynomial algorithm |21| . It is 
worth doing this step, since it might save some compntational time and space, bnt it is 
not necessary for the rnnning of the algorithm. 

Based on the above, each reaction graph can be nniqnely determined if it is known 
which non-core reactions of the dense realization take place in the reaction network. Thns, 
we represent the reaction graphs by binary seqnences of length N = \E{Gd) \ Fc|. 
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In order to define the binary sequences we fix an ordering of the non-core edges. Let 
Cj denote the ith edge. If i? is a binary sequence, then let R[i] denote the ith coordinate 
of the sequence and Gr be the reaction graph described by R. If a realization can be 
decoded by the binary sequence i?, then 

fe e 

e e E{Gr) < or (19) 

[d* e {1,... ,iV} e = ei,R[i] = l 

From now on the term ‘sequence’ will refer to such a binary sequence of length N . 
The sequence representing the dense realization (with all coordinates equal to 1) will be 
denoted by D. 

For the efficient operation of the algorithm, appropriate data structures are needed as 
well. The discovered graph structures are stored in a binary array of size 2^ called Exist, 
where the indices of the helds are the sequences as binary numbers. At the beginning the 
values in every held are zero, and after the computation the value of held Exist[R\ is 1 if 
and only if there is a linearly conjugate realization described by the sequence R. 

We also need N + 1 stacks, indexed from 0 to N. The kth stack is referred to as S{k). 
During the computation, sequences are temporarily stored in these stacks, following the 
rule: the sequence R describing a linearly conjugate realization might be in stack S{k) 
if and only if there are exactly k coordinates of R which are equal to 1, i.e. exactly k 
reactions take place in the realization. 

At the beginning all stacks are empty, but during the running of the algorithm we 
push in and pop out sequences from them. The command ‘push R into *S'(fc)’ pushes the 
sequence R into the stack S{k), and the command ‘pop S{ky pops a sequence out from 
S{k) and returns it. (It makes no diherence, in what order are the elements of the stacks 
popped out, but by the dehnition of the data structure the sequence pushed in last will 
be popped out hrst.) The number of sequences in stack S{k) are denoted by size.S'(fc), 
and the number of coordinates equal to 1 in the sequence R are referred to as e{R). 

Within the algorithm the following procedure is used repeatedly: 

FindLinConjWithoutEdge(M, y, i?, i) computes a constrained dense linearly conju¬ 
gate realization of the kinetic system with coefficient matrix M and complex composition 
matrix Y . The additional inputs R and i are a sequence encoding the input reaction 
graph structure, and an integer index, respectively. The procedure returns a sequence U 
encoding the graph structure of the computed realization such that Gr is a subgraph of 
Gr and U[i\ = 0. If there is no such realization, then —1 is returned. This computation 
can be carried out in polynomial time as it is described in |3n) . 
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Algorithm 1 Determines all reaction graphs describing linearly conjugate realizations 

1: procedure Linearly conjugate graph structures(M, F) 

2: push D into S{N) 

3: Exist[D]:=l 

4: for A; = iV to 1 do 

5: while size.S{k) > 0 do 

6: R:= pop S{k) 

7 : for i = 1 to iV do 

8: if R[i] = 1 then 

9: U := FindLinConjWithoutEdge(M, F, i?, i) 

10: if f/ > 0 and Exist\U] = 0 then 

11: Exist[U]:= 1 

12: push U into S{e{U)) 

13: end if 

14: end if 

15: end for 

16: Print R 

17: end while 

18: end for 

19: end procedure 


Proposition 4.1. For any kinetic system and any suitable fixed set of complexes all the 
possible reaction graphs describing linearly conjugate realizations can be computed after 
finitely many steps by Algorithm 1. The whole computation might last until exponential 
time depending on the number of different reaction graphs, but the time elapsed between 
the displaying of two linearly conjugate realizations is always polynomial. 

Proof. Let us assume that there is a sequence W which is not returned by the algorithm, 
but it describes a linearly conjugate realization of the kinetic system. Let R be another 
realization, which was computed by the algorithm and Gw is a subgraph of Gr, i.e. 
R[i] = 1 if W[i] = 1 for alH G {1,..., iV}. The sequence D fulfils this property for each 
IF, but if there is more than one such sequence, then let us chose R to be the one with 
minimum number of coordinates equal to 1. It follows from the definition of the sequences 
IF and R that there must be an index j so that W[j] = 0 and R[j] = 1. (If there is more 
than one such index, then let us chose the smallest one.) During the computation there 
is a step (line 9) when we apply procedure FindLinConjWithoutEdge(M,F,i?,j) to 
compute the sequence U, which is a dense realization with the properties: U[j] = 0 and 
Gu is a subgraph of Gr. According to the properties of IF, it fulfils the constraints of 
the optimization problem specified in the procedure, therefore it cannot be infeasible and 
Gw rnust be a subgraph of Gr. But it leads to a contradiction, since if IF is equal to U, 
then IF is returned by the algorithm, and if they are not equal, then R is not minimal. 

A sequence R popped out from stack S{k) is printed out after all its coordinates 
are examined. This computation requires the application of the procedure FindLinCon- 
jWithoutEdge(M,F,i?,i) k times, with some additional minor computation. From the 
properties of the procedure it follows that the computation between displaying two con¬ 
secutive realizations can be performed in polynomial time. 

In each stack there might be finitely many sequences (in stack S{k) at most (^)) and 
in case of each sequence, only finitely many calls of FindLinConjWithoutEdge have 
to be performed, therefore the whole computation can be performed in finite time. □ 
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It is important to see that we compute the entries of the corresponding matrices Ak 
and T~^ in each step of Algorithm 1 according to Eqs. (12)-(16), but do not store them. 
This is possible because to set up the constraints for the subsequent computations, it 
is enough to know only the structure (i.e. the zero and non-zero entries) of matrix 
corresponding to the realization that is popped out of the actual stack. If one wants to 
store and/or display the actual reaction rate coefficients and transformation parameters, 
it is possible by slightly modifying the procedure FindLinConjWithoutEdge to return 
not only the resulting reaction graph encoded by the binary sequence [/, but also the 
matrices T and A^. 


4.1 Parallelization of the algorithm 


Due to the possible large number of different reaction graphs even in the case of rel¬ 
atively small networks (see Example 2 in Section [^, it is worth studying the parallel 
implementation of the proposed algorithm. This may dramatically improve the computa¬ 
tional performance as it has been shown in the case of several other fundamental problems 
(see, e.g. [33l [3l]). 

For the sequences R in each stack S{k) and for all the indices i E {1,..., N} the results 
coming from the executions of the procedure FindLinConjWithoutEdge(M, Y, R, i) do 
not have any effect on each other, consequently these might be computed in a parallel 
way. The obtained sequences are pushed into stacks with indices smaller than k, but only 
the ones that have not been found earlier. It follows that the procedure might be applied 
in parallel for sequences from different stacks as well, since the repetition of the same 
computation is avoided by the application of the array Exist. 

In case of dynamically equivalent realizations it is possible to do even more paralleliza¬ 
tion. These are special linearly conjugate realizations, where the transformation matrix 


T is the unit matrix, the variables are the entries of matrix Ak, and Equation (14) deter¬ 
mining dynamical equivalence can be written in a more simple form: 


Y ■Ak = M (20) 

It is easy to see that the values of the variables in the jth column of matrix Ak depend 
only on the parameters in the jth column of matrix M and on the entries of matrix 
Y . Therefore the columns of Ak might be computed parallely, and any dynamically 
equivalent realization can be determined by choosing a possible solution in case of each 
column and build the Kirchhoff matrix of the realization from them. It follows that the 
number of dynamically equivalent realizations, that describe different reaction graphs, is 
the product of the numbers of the possible columns. 

The super-structure property of dense realizations is inherited by the columns of the 
matrix Ak, therefore we can use the same algorithm for these as we used for linearly 
conjugate realizations. 

In this algorithm it is better to determine the ordering of non-core edges according 
to columns. Let us denote the number of non-core edges in column j by Nj, and the 
sequence describing the jth column of the dense dynamically equivalent realization by 
Dj. The stacks are also needed to be dehned separately for each column. The sequence 
Rj representing a jth column gets stored in stack Sj{k) if and only if the number of 
coordinates equal to 1, denoted by e{Rj), is exactly k. 

We also need a two-dimensional binary array denoted by ExistColumn[j, Rj] to store 
the computed sequences in case of each column. The hrst index refers to the column and 


11 



the second index is the sequence as a binary number. At the beginning all coordinates are 
equal to zero. 

The applied procedures are as follows: 

• DyneqColumnWithoutEdge(M, Y, j, Rj, i) computes the jth column of the Kirch- 
hoff matrix describing a constrained dense dynamically equivalent realization of the 
kinetic system with coefficient matrix M and complex composition matrix Y. The 
constraints are determined by the two last inputs, a sequence Rj and an integer 
index i. The procedure returns a sequence Uj representing a jth column so that 
Uj[l] = 0 if Rj[l] = 0 for all Z G {1,... A^} and Uj[i] = 0. If there is no such column, 
then —1 is returned. This computation can be performed in polynomial time. 

• Bui\dAk{ExistColumn) builds all possible dynamically equivalent realizations from 
the sequence parts in ExistColumn and saves them in the array Exist. 


Algorithm 2 Determines all reaction graphs describing dynamically equivalent realiza¬ 
tions applying parallelization 

1 

procedure Dynamically equivalent graph structures(M, Y) 


2 

for j = 1 to m do 


3 

push Dj into Sj{Nj) 


4 

for k = Nj to 1 do 


5 

while size.Sj{k) > 0 do 


6 

Rj.= pop Sj{k) 


7 

for i = 1 to Nj do 


8 

if Rj[i] = 1 then 


9 

Uj := DyneqColumnWithoutEdge(M, Y, j, i) 


10 

if Uj > 0 and ExistColumn[j, Uj] = 0 then 


11 

ExistColumn[i, Uj\-.=1 


12 

push Uj into Sj{e{Uj)) 


13 

end if 


14 

end if 


15 

end for 


16 

end while 


17 

end for 


18 

end for 


19 

BuildAfc [ExistC olumn) 


20 

end procedure 



5 Examples 


In this section we demonstrate the operation of our algorithm on kinetic systems with 


a small (Section 5.1) and a slightly bigger (Section 5.2) set of complexes. The algorithm 


was implemented in MATLAB pa using the YALMIP modelling language |36) . 

We will see that the number of possible reaction graphs describing linearly conjugate 
realizations grows very fast depending on the number of complexes. 
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5.1 Example 1 

In this example taken from ra. we examine the kinetic system described by the following 
dynamical equations: 


Xi = 3^1 ■ x\ — k2- x\ 

X2 = —3ki ■ xl + k2 ■ xl 

According to the monomials there are (at least) two species, S = {Xi, X 2 }, and we £x the 
set of complexes to be C = {Ci, C 2 , Ca}, where Ci = 8 X 2 , C 2 = 3Xi, and C 3 = 2Xi + X 2 . 
Based on the above, the inputs of the algorithm - the matrices Y and M - are as follows: 

3fci —^2 0 
—3/ci ^2 0 

For the numerical computations, the parameter values ki = 1 and ^2 = 2 were used. 
As the result of the algorithm we get 18 different sequences/reaction graphs. This small 
example is special in the sense that the sets of different reaction graphs corresponding 
to dynamically equivalent and linearly conjugate realizations are the same, since the 
computed transformation matrix T was the unit matrix in each case. Using the numerical 
results, it was easy to symbolically solve the equations for dynamical equivalence, therefore 
we can give the computed reaction rate coefficients as functions of ki and ^ 2 - 

The reaction graphs are denoted by Gi,, Gig and are presented with suitable re¬ 
action rate coefficients in Figure [Tj From the computation it follows that there are two 
reaction rate coefficients k^i and ^32 which do not depend on the input parameters, just 
on each other, and the reactions determined by these might together be present or 
non-present in the reaction network. Therefore, a nonnegative parameter p is applied to 
determine the values of these coefficients. We get the reaction graphs Gi,..., Gg if the 
parameter p is positive, and if it is zero then we get the reaction graphs Gio,..., Gig. 
The reaction graph Gi (the complete directed graph) describes the dense realization, 
and consequently all other reaction graphs are subgraphs of it (not considering the edge 
weights). 


Y = 


0 3 2 
3 0 1 


M = 
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5.2 Example 2 


The purpose of this example is to show the possible large number of structurally different 
linearly conjugate realizations even in the case of a relatively small kinetic system. The 
reaction network examined in this section was published in |3H] as example Al. In the 
original article it is given by the following realization, described by the Kirchhoff matrix 
Ak and reaction graph shown in Figure 


Cl 



Ak = 


-ki 

0 

ki 

0 

0 

0 


k2 

-k2 

0 

0 

0 

0 


0 

^3 

-^3 

0 

0 

0 


0 

0 

^4 

-^4 

0 

0 


0 0 
0 0 
0 0 
0 0 


-h 

h 


Figure 2: The reaction graph of Example 2 

In the reaction network there are two species, S = {Xi,X 2 } and six complexes, C = 
\C\ = 0,(^2 = Xi^C^ = X 2 ,C/^ = 2Xi,C^ = 2Xi + X 2 iCq = According to the 

dehnitions, the complex composition matrix Y and the matrix M = Y ■ Ak oi coefficients 
are as follows: 


Y = 


0 

1 

0 

2 

2 

3' 

M = 

'0 

-k2 ks 

—2ki 

fcs 0 

0 

0 

1 

0 

1 

0 

ki 

0 -ks 

ki 

-fcs 0_ 


The reaction rate coefficients used in the computations were the same as in j3H], namely: 
ki = 1, k 2 = 1, ^3 = 0.05, k^ = 0.1, fcs = 0.1. With these parameter values the system 
shows oscillatory behaviour. In the dense linearly conjugate realization shown in Figure 
there are 19 reactions, and it can be given by the matrices and T‘^ as: 


Cl 



Ad _ 
— 


C3 


80 

1.167 

10^ 

3.083 

3.333 

106 

0.333 

O' 

0 

-2 

10^ 

0.5 

5 

10^ 

0.5 

0 

80 


0 

-4.25 


4 

0.5 

0 

0 

5 

10® 

0.25 

-2 

10" 

1 

0 

0 


0 

0.25 


4 

-8.5 

0 

0 

3.333 

10® 

0.167 

1.167 

10" 

6.167 

0 


Figure 3: The reaction graph of the dense 
linearly conjugate realization 


(yd)-i ^ 


40 0 

0 80 


Our algorithm returned as many as 17160 different reaction graphs describing linearly 
conjugate realizations of this kinetic system, all of which can be found in the electronic 
supplement available at: 
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http://daedalus.scl.sztaki.hu/PCRG/works/publications/Ex2_AllRealSuppl.pdf 

Out of these, 17154 can be described by a weakly connected reaction graph, while 6 have 
disconnected reaction graphs, with the same linkage classes. Since this property can be 
ensured by linear constraints (the edges between the linkage classes are excluded), accord¬ 
ing to Proposition |2.8| the realization having the maximum number of edges determines a 
super-structure among realizations obeying the same constraints. This constrained dense 
realization can be described by the matrices and T^'^, and its reaction graph is shown 
in Figure 



-50 

1.25 

w 

0.625 

0 

0 

0‘ 

0 

-2.5 

10" 

1.25 

0 

0 

0 

50 


0 

-2.5 

5 

0 

0 

0 

1.25 

10" 

0.625 

-5 

0 

0 

0 


0 

0 

0 

-5 

0 

0 


0 

0 

0 

5 

0 



1 _ 

50 O' 
0 50 





Figure 4: The reaction graph of the dense 
linearly conjugate realization having two 
linkage classes 

It also turned out from the computations that in this case the sparse realization is 
unique, and it is the initial network shown in Figure]^ The distribution of the computed 
reaction graphs over the number of reactions is shown in Fig. 



Number of reactions contained in the realizations 


Figure 5: Number of different reaction graphs with given numbers of reactions (directed 
edges) in the case of Example 2 
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Figures and show the solutions of the original and the dense linearly conjugate 
realizations from consistent initial conditions. The oscillatory behaviour is clearly shown, 
and one can check from the results that x{t) = (T^)“^ • x{t) holds for alH > 0. 



Figure 6: Solution of the original system in Example 2 defined by the matrix pair [Y, Ak) 
from the initial conditions a;(0) = [1 2]^ 



Figure 7; Solution of the linearly conjugate system in Example 2 dehned by the matrix 
pair {Y,A'l) from the initial conditions a:(0) = ■ a;(0) = [40 160]^ 

5.3 Computation results of the parallel implementation 

We tested the parallel implementation of Algorithm 1 on a workstation with two 2.60GHz 
Xeon (E5-2650 v2) processors with 32 Gb RAM (DDRS 1600 MHz, 0.6ns). The imple¬ 
mentation was written in Matlab R2013a using the built-in parallelization toolbox. 

Let L denote the number of threads. The tests were carried out with L G {1, 2, 4, 6, 
8, 10, 12}, where these threads are working on the innermost loop (lines 7 to 15) of the 
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algorithm. We can easily conclude that the most time-consuming step is the calling of 
procedure FindLinConjWithoutEdge which required on average 1.9 ms computation 
time with 0.6 ms standard deviation in case of Example 2 (calculated considering 130 000 
executions). Figure shows the total execution times in case of both examples with 
different numbers of threads on log-log scales. 



Figure 8: Overall execution times of Examples 1 and 2 with different numbers of threads 

One can expect that the computation of N different edge exclusions by using the 
procedure FindLinConjWithoutEdge reaches maximum efficiency when L = N holds. 
It can be seen that in the case of Example 1 (where iV = 6) we reach the maximum 
efficiency indeed, and for a larger number of threads we get slightly longer computation 
times. The improving efficiency with the larger number of threads can be seen in the case 
of Example 2, where the value of N was 19. The results show in the studied thread-range 
that by doubling the number of threads, the total execution time approximately gets 
halved. 

In case of Example 2 we have recorded the computation times separately for the 
different stacks as well. Here the computation time means the time elapsed from popping 
out the first element from a given stack until the end of processing the last element in 
the same stack. The obtained results are shown in Table Since the sparse realization 
contains five reactions, the procedure FindLinConjWithoutEdge was not called for 
sequences with fewer than five nonzero coordinates. This is why stacks with index smaller 
than 5 are not included in the table. 


6 Conclusions 

An algorithm was proposed in this paper for computing all structurally different reaction 
graphs describing linearly conjugate realizations of a given kinetic polynomial system. 
To the best of the authors’ knowledge this method is the first provably correct solution 
for the exhaustive search of CRN structures realizing a given dynamics. The inputs of the 
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Number of threads 



1 

2 

4 

6 

8 

10 

12 


5(5) 

17.667 

5.7767 

2.6061 

1.9417 

1.0794 

0.91359 

0.83738 


5(6) 

149.11 

42.734 

19.173 

16.935 

8.763 

7.3073 

6.5739 


5(7) 

762.26 

234.79 

87.411 

100.92 

42.237 

35.511 

31.838 


5(8) 

4458.7 

1233.5 

482.81 

532.03 

236.23 

197.91 

177.8 


5(9) 

20355 

6247.5 

2503.8 

2356.1 

1278.1 

891.4 

794.43 


5(10) 

59492 

17246 

6921.6 

6769.3 

3511 

2536.5 

2248.6 


5(11) 

1.0534e+05 

33706 

13570 

11661 

5982.5 

5338.1 

3878.8 

Stacks 

5(12) 

1.1675e+05 

36972 

20753 

11958 

6507.2 

5825.6 

4253.1 


5(13) 

83201 

25604 

14972 

8653.4 

4572.2 

4085.2 

3770.1 


5(14) 

35556 

10750 

5960.5 

3755.5 

2126.5 

1958.5 

1858.1 


5(15) 

9283.8 

3501.6 

1737.8 

1197.8 

787.78 

749.61 

723.94 


5(16) 

2015 

898.08 

478.54 

353.27 

259.49 

246.49 

234.34 


5(17) 

397.5 

205.36 

125.82 

81.447 

75.907 

57.191 

55.652 


5(18) 

57.236 

29.231 

17.953 

11.173 

11.078 

8.6908 

8.0421 


5(19) 

4.9291 

3.3659 

2.7989 

2.0673 

1.9249 

1.8916 

1.7575 


Table 1; Processing times (in seconds) of the individual stacks in case of Example 2 for 
different numbers of threads 


algorithm are the complex composition matrix and the coefficient matrix of the studied 
polynomial system. The output is the set of all possible reaction graphs encoded by binary 
sequences. The correctness of the method is proved using a recent result saying that the 
linearly constrained dense realization determines a super-structure among all realizations 
fulhlling the same constraints [30]. Although exponentially many different reaction graphs 
may exist, it is shown that polynomial time is elapsed between displaying (storing) two 
consecutive reaction graphs. The computation starts with the determination of the dense 
realization and different stacks are maintained for storing realizations with the same num¬ 
ber of reactions. The number of stacks depends linearly on the number of reactions in the 
dense realization, although the entire bookkeeping may require exponential storage space 
due to the possible large number of different structures. This organization allows that the 
optimization tasks for processing the actual realization (i.e. computing its constrained 
immediate ‘successors’) can be implemented parallely. The applicability of the algorithm 
is illustrated on two examples taken from the literature. The numerical results show that 
parallel implementation indeed improves efficacy of the computation. 
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